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Abstract - In this letter we consider the phase diffusion of a harmonically driven undamped pen¬ 
dulum and show that it is anomalous in the strong sense. The role played by the fractal properties 
of the phase space is highlighted, providing an illustration of the link between deterministic chaos 
and anomalous transport. Finally, we build a stochastic model which reproduces most properties 
of the original Hamiltonian system by alternating ballistic flights and random diffusion. 


Anomalous transport is a phenomenon relevant to 
a wide range of complex systems and has recently at¬ 
tracted the attention of researchers with a mathematical, 
physical, chemical, biological and socio-economical back¬ 
ground. The list of diffusion processes where violations 
of the hypotheses of the central limit theorem lead to de¬ 
partures from the normal long-time asymptotic behaviour 
{x^{t)) ~ t, where x(t) is a generic dynamical variable, is 
huge and ever-growing due to an intense theoretical and 
experimental research effort [^. Here we focus on strong 
anomalous diffusion [^, defined by the property 

( 1 ) 


where the function v{q) is not constant. Such behaviour 
has been detected in numerical studies over a variety of 
systems: without any claim of completeness, we can cite 
intermittent ID maps , running sandpile models ID , in¬ 
finite horizon and polygonal billiards, cold atoms in 
optical lattices [S] and ID inhomogeneous materials [^. 
The emergence of strong anomalous diffusion has been 
analytically investigated in the context of both stochas¬ 
tic [lOpT and Hamiltonian models. The first (and, at 
the moment, only) experimental measure of strong anoma¬ 
lous diffusion has been very recently obtained by tracking 
polymeric particles in living cells 13 . In this letter we 


will investigate the diffusion properties of the phase of a 
harmonically driven pendulum. This system is representa¬ 
tive of the features of generic Hamiltonian continuous-time 
dynamics but at the same time it is simple enough to be 
considered a paradigmatic model for anomalous transport. 


Chaotic pendulum. — 

familiar nonlinear system 


The pendulum is the most 
14 and its study goes back to 


the very first appearance of modern (as opposed to an¬ 
cient) physics. When subject to a suitable periodic driv¬ 
ing force, the pendulum displays chaotic behaviour 
one can then investigate its phase diffusion properties 


15 


16 


by choosing many initial conditions in a tiny region of 
the phase space and observing how the deterministic tra¬ 
jectories spread out. The periodically driven undamped 
pendulum equation we consider here 


( 2 ) 


9 + sin 0 = 7 sin(a;t) 


has been studied previously by the authors of refs. |17f|20| , 
who report instances of both regular (in the case 7 = 1.2, 
uj = 0.1) and anomalous (for 7 = 1.2, uj = 0.8) diffusion. 
Actually, we partially correct here those statements, show¬ 
ing explicitly that the diffusion is strongly anomalous in 
both cases. 

There are good reasons to believe that strong anomalous 
diffusion is rather the rule than the exception in Hamilto¬ 
nian systems characterised by a mixed phase space, where 


regular islands are surrounded by the chaotic sea 21-23 
A chaotic trajectory cannot enter an island, which is 
formed by KAM-tori, yet it can temporarily behave like 
a regular motion due to the presence of cantori 24 all 
around the island. Cantori are fractal objects that look 


like closed curves with an infinite number of gaps 25 ; they 


act on trajectories as quasi-traps, i.e., they cause the ve¬ 
locity to remain constant during sizeable periods of time, 
whereas a purely chaotic motion would lead to nearly un- 
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Fig. 1: Function /iq('r) as defined by eq. 0 for different values of q. In order to show that the log-periodic oscillations are not 
an artifact due to poor sampling, we have computed 10® additional (shorter) traces, which contribute to the thick lines. 


correlated displacements. The kinetics is then similar to 
that of a Levy walk 26,27, a process where the system 


randomly select a velocity and maintains it for a random 
time before renewing the choice. Such process is known 
to give rise to strong anomalous diffusion [^. The pre¬ 
vious quite general argument suggests that there is noth¬ 
ing special in our choice of the parameters ( 7 , w) and that 
the same conclusions would be obtained by inspecting any 
other value compatible with chaotic behaviour. 

In the present work we numerically integrate eq. 0 
using a Runge-Kutta 4th order algorithm [^. The ad¬ 
equacy of the chosen integration step {5t = 0.005 for 
w = 0.1 and 5t = 0.001 for oj = 0.8) was confirmed by 
checking that no statistical property of the trajectories is 
signihcantly altered when the integration step is halved. 
For each value of w, about 5 x 10'^ initial conditions (0Oj ^0) 
are randomly chosen in the region —0.05 < 0n < 0.05, 
—0.05 < 00 < 0.05. The time t appearing in eq. (0 is mea¬ 
sured in units of the proper frequency of the pendulum, 
but in order to compare behaviours for different values 
of ui it is useful to present results in terms of the nor¬ 
malised time r, measured in cycles of the external force: 
T = ujtl{27r). We follow the evolution of our sample tra¬ 
jectories along 10 ® such cycles. 

Phase diffusion properties. — There are a few stan¬ 
dard ways to characterise anomalous diffusion 


29 . One 


is to evaluate the asymptotic exponent qv{q) of fractional 
moments mq{T): 


mq(T) = (| 0 (t)|®) ^ for 


(3) 


The study of rnqir) over several orders of magnitude in 
time reveals small but not negligible departures from a 
straight power law (see hg. , which are more clearly vi¬ 
sualised if we consider the quantity 


= r—Inm,(r) , 


(4) 



Fig. 2: Statistical momenta m, as a function of the normalised 
time r, for ui = 0.8. (The case a; = 0.1 yields very similar 
results.) 


showed in fig. Similar log-periodic oscillations have 
been related to the presence of self-similar structures in 
the phase space and occur frequently in the litera¬ 
ture about diffusion on fractals 31-33 . For our purposes 
here, the existence of these oscillations implies that the 
asymptotic exponent qi'{q) markedly depends on the time 
interval over which it is estimated. 


With this caveat in mind, we select the largest time 
window permitted by our data, r € [ 10 ^, 10 ®], and use 
that region to evaluate qv{q). The results for our chaotic 
pendulum are reported in fig. If we fit a power law 
to mqij) over another time interval, the precise values of 
the exponents change, but the fact that oscillations for 
different q's are not in phase tells us that for all but a 
carefully picked set of choices we find the same qualita¬ 
tive behaviour, where we easily recognise the hallmarks 
of strong anomalous diffusion: moments with low index 
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Fig. 3: Asymptotic exponent of the fractional moments as defined by eq. HI- Strong anomalous diffusion is made evident for 
both values of uj by the existence of two regimes: superdiffusive and ballistic for lower and higher values of q, respectively. 
Normal diffusion would correspond to the straight line through the origin with slope 1/2. 


q are dominated by the most likely trajectories and ex¬ 
hibit superdiffusion due to the aforementioned mixed na¬ 
ture of the phase space; high-g moments, on the other 
hand, are dominated by the ballistic behaviour of the ex¬ 
tremal traces, those that have spent most time around 
the regular regions and have therefore traveled furthest (a 
detailed account of the behaviour of the solutions 0{t) is 
given in the next section). 

Another way to visualise the nature of the diffusion is 
to consider the probability density function P{9,t) for the 
phase 9 at time r. In the case of a normal diffusion process 
(t/(q) = 1/2 for any q), we would expect to be able to find 
an exponent rj and a function $ such that 

= (5) 

with rj = 1/2. For weak anomalous diffusion {^{q) = v 
constant but different from 1 / 2 ) the scaling form above 
still holds, but in general r] ^ v. In the case of strong 
anomalous diffusion {iy{q) not constant), however, it is not 
possible to find 77 such that eq. ([^ is satisfied for all the 
values of 9: the collapse of the curves P{9, r) for different 
r is limited to the central part of the distribution and 
breaks down in the large -0 regime (see fig. dominated 
by ballistic events for which the scaling variable would 
rather be 9/t. 

In ref. it is showed that for a class of continuous¬ 
time random walks it is 77 = limg_,,o v^q)- Here we compare 
the slope 0.66 from fig. with 77 = 0.61 which gives the 
best superposition of the curves <1>(0 /t'') for w = 0.8. In 
the case w = 0.1 the same comparison can be drawn be¬ 
tween 0.59 from the right panel of fig.|^and the empirically 
determined 77 = 0.54. 

Phase space and trajectories. — In the previous 
comments about the diffusion properties we have repeat¬ 
edly emphasised the role played by the phase space ( 0 , 0 ); 



O/r” 

Fig. 4: Rescaled probability distribution >1>(0 /t’’) for cj = 0.8. 
The value rj = 0.61 was found by empirically optimising the 
collapse and is quite close to lim^^o ~ 0 . 66 . 


let us now try to visualise its structure and how such struc¬ 
ture is reflected by the trajectories 0 (t). In the left panel 
of fig. [^a solution of eq. Q obtained for oj = 0.8 and ini¬ 
tial condition close to (0, 0) is represented by its Poincare 
section: the phase 0 and angular velocity 0 are recorded at 
the end of every cycle of the external force. On the right 
panel of the same figure we draw a segment of the 0 (r) 
graph of the same solution. We clearly recognise the typ¬ 
ical features of a mixed phase space: even though regular 
regions are inaccessible to a chaotic trajectory (and are 
therefore indicated by blank spaces in the Poincare sec¬ 
tion), their presence influence the dynamics by creating 
traps whence, due to barriers formed by cantori, it takes 
a long time for the system to get out. While trapped in 
these regions the pendulum maintains a constant veloc¬ 
ity and as a consequence large variations of 0 occur, in 
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Fig. 5: On the left panel, the Poincare section of one trajectory, obtained with ca = 0.8; on the right panel, a segment of the 
same trajectory. The colours highlight regions of the phase space where the dynamics gets trapped, so that the velocity stays 
approximately constant for a macroscopically long time. 


contrast with the random walk that takes place while the 
system is exploring the so called chaotic sea. 


The deterministic trajectory in fig. right panel, re¬ 
sembles the result of a stochastic process belonging to 
the family of Levy walks 34 . Following the method of 
ref. 27 , we can explicitly construct such process. As a 


first step, we need to identify a small number of sticky 
regions, which is easily accomplished by considering the 
distribution of A9 = 9{t + At) — 6{t) for a suitable fixed 
value of At. In the case uj = 0.8 the choice At = 50 per¬ 
mits the easy identification of six peaks symmetrically ar¬ 
ranged around 0 (see hg. [^, which match the six coloured 
regions in fig. We can then perceive the excursions in 
the chaotic sea, during which 9 undergoes normal diffusion 
with zero mean, as pauses between ballistic flights, which 
happen when the system is trapped into one of the six 
regions [i = 1,2,3). See the caption of fig. for the 
correspondence between colour codes and typical mean ve¬ 
locities A9/ At (not to be confused with the instantaneous 
velocities 9 = d9/dt). 


The time between flights is exponentially distributed 
ipdt) = T“^e“*/'^‘= (in our case « 270 cycles of the ex¬ 
ternal force), while the duration of each flight follows a 
power-law tpid) ^ To each ballistic region we can 

therefore associate three quantities: 1) a winding number 
Wi = A9/ (2ttAt), which is the velocity maintained during 
the flight, 2) the probability Pi of entering region i from 
the chaotic sea and 3) the exponent ai that characterises 
the residence time within the region. We have already 
identified the winding numbers as the peaks in hg. as 
for the probabilities pi and residence time exponents Ofi, 
we divide the time in intervals of duration At = 10 cycles 
of the external force and associate each interval to one bal¬ 
listic region or to the chaotic sea according to the mean 
velocity falling into one of the peaks or not. The parame¬ 
ters that better describe the phase diffusion of our chaotic 



Fig. 6: Distribution of A9/ At obtained for cj = 0.8 with At = 
50. The dashed line is a Gaussian approximation to the central 
part. The typical velocities of the six coloured regions in fig.[^ 
are signaled by the peaks: from left to right (green), B 2 
(magenta), B^ (cyan), B^ (blue), 5^ (yellow) and B^ (red). 


pendulum for w = 0.8 are summarised in tablej^ Observe 
that for the regions we report two exponents, which 
dominate the behaviour of ip 2 it) for t < 80 ( 0:2 = 0.5) and 
for t > 80 [a '2 = 2.5). Of course, the identification of just 
six ballistic regions is a simplification of the real complex¬ 
ity of the phase space structure: it may well be that the 
velocity-based approach we are adopting here simply does 
not have enough resolving power to distinguish regions 
with similar winding number but quite different residence 
time distributions, and this is reflected in the lack of a 
definite power law for regions B^. 

Note that previous studies of the phase space structure 
of this same system 


regions Bf 


20 


focused on the extreme velocity 
overlooking the role played by B^ 
which we on the contrary find determinant to describe the 


and B 


3 ’ 
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Table 1: Parameters of the kinetic model that mimics the dif¬ 
fusion properties of the phase of the chaotic pendulum for 
cj = 0.8. The system diffuses during an exponentially dis¬ 
tributed time, and then enters with probability p; one of the 
six ballistic regions characterised by winding number Wi and 
residence time exponent q^. 


the figure is not calculated from a fit, whose result would 
depend on how many points are taken into account, but 
from tracing the straight line between the first two data 
available: q = {) and q — 0.1. This gives another check 
that the kinetic model we have built is consistent with the 
results of our numerical analysis. 



Wi 

Pr 

Oii 


±1.0 

0.296 

1.5 

B^ 

±4.0 

0.197 

0.5 / 2.5 

Bi 

±4.2 

0.007 

0.5 


superdiffusive regime. Actually, from the data in table 
we can estimate a theoretical value for the diffusion expo¬ 
nent 2v{2). As the presence of exponentially distributed 
waiting times between ballistic flights has no effect on the 
mean squared displacement 35 , we can use the scaling 


laws 36 37 




t^ 

0 

< 

a < 


a 

= 

1 

^3-a 

1 

< 

a < 

tint 

a 

= 

2 

t 

a 

> 

2 


( 6 ) 


Conclusions. — We have presented numerical evi¬ 
dence that the phase of the harmonically driven undamped 
pendulum of eq. ([^ undergoes strong anomalous diffusion 
in both the cases we have analysed (7 = 1.2, a; = 0.1 or 
0.8). We have attributed such phenomenon to the mixed 
nature of the phase space, with sticky regions near the 
regular solutions that act as dynamical traps, so that the 
time evolution of 0 resembles that of a stochastic process 
where standard diffusion is interrupted by ballistic flights 
with power-law distributed duration. This slow decay of 
the flight-time distribution is responsible for the superdif¬ 
fusive behaviour, while the existence of extremal flight ve¬ 
locities determines the strongly anomalous character of the 
diffusion process, as can be seen in analytical calculations 
based on a Levy walk model [H 34 . 


With the idea of suggesting a system where strong 
anomalous diffusion could be experimentally observed, we 
plan in the future to investigate to what extent these fea¬ 
tures survive the introduction of a friction term. 


derived for a Levy walk characterised by the flight time 
distribution 4’{t) ~ 

The probability distribution of the residence time in B 2 
falls off rapidly (exponent 2.5) for large t, so we can assume 
these regions not to play a dominant role in determining 
the diffusion exponent. On the other hand, even if regions 
host the longest ballistic excursions, they are rarely 
accessed: we expect them to influence the extreme trajec¬ 
tories (relevant to moments of larger order), but the mean 
squared displacement should be dominated by the more 
frequently visited B^ regions. If we adopt ai = 1.5 as 
the exponent that better characterises the phase diffusion 
process, then eq. § yields {OitY) ~ which is in fact 
the value we observe (see fig.[^. 

The same analysis can be applied to the case a; = 0.1 
with similar results: the main difference is that the kinetic 
model is simpler as we need to identify only two ballistic 
regions B^, characterised by winding number w = ±140 
rounds per cycle of the external force and by the exponent 
a = 1.7. Once again the expectation for the diffusion 
exponent based on eq. matches the value 1.3 from our 
numerical test. 

Another quantity that can be analytically computed in 
a Levy walk model is the behaviour of v{q) for small g, 

Mmg^o iy{q) = 1/a. 


for which we have the relation 10 


Taking ai = 1.5 as the dominant exponent for w = 0.8 and 
a = 1.7 for uj = 0.1 we would expect therefore that i/{q) = 
0.67 and 0.59, respectively, when q is small, in remarkable 
agreement with fig. Note that the slope for small q in 


Part of this research took place at the Galileo Galilei 
Institute for Theoretical Physics in Arcetri, during 
the INFN-funded summer 2014 workshop Advances in 
Nonequilibrium Statistical Mechanics. 
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